%% Finite Slit Diffraction

lambda=1;
k=2*pi/lambda
a=20;
w=10;

Xmax=200
Ymax=50

%Use 'N' sources to represent the slit
N=10;

Svec=linspace(-a/2,a/2,N)
WaveAmp=0;
for ss=1:N
    [X,Y]=meshgrid((0:.1:Xmax),(-Ymax+Svec(ss):.1:Ymax+Svec(ss)));
    R=sqrt(X.^2+Y.^2);
    Wave=exp(i*(k*R));% ./R1; %From Huygens Principle
    WaveAmp=WaveAmp+Wave;
end

figure
imshow(WaveAmp,[])

%Now you could take the time average like we did on April 10



